clear;
close all;
clc;

% Output Selection

report       = 'T5';          % 'T5', 'T6', 'T7' 
markup       = '1.15';        % '1.05', '1.15', '1.25', '1.35'
experiment   = 'efficient';   % 'efficient', 'uniform', 'sizedep', 'entry'

% Options

set(groot, 'DefaultAxesLineWidth', 1.5);
set(groot, 'DefaultLineLineWidth', 3);
set(groot, 'DefaultAxesTickLabelInterpreter','latex'); 
set(groot, 'DefaultLegendInterpreter','latex');
set(groot, 'DefaultAxesFontSize', 22);
   
options   = optimoptions('fsolve', 'FiniteDifferenceStepSize', 1e-4, 'Display', 'none', 'TolX', 1e-11, 'TolFun', 1e-11, 'UseParallel', false);

% Calibrated Parameters

switch markup 

    case '1.05'

        Markuptarget = 1.05;                      % our target for aggregate markup
        
        N            = 414.6473;                  % mean number of firms in sector
        p.sigma      = 1.6220;
        p.gamma      = 59.6930; 
        p.xi         = 28.0803; 

    case '1.15'

        Markuptarget = 1.15;                      % our target for aggregate markup
        
        N            = 358.61;                    % mean number of firms in sector
        p.sigma      = 1.3506;
        p.gamma      = 12.758; 
        p.xi         = 8.5109; 

    case '1.25'    

        Markuptarget = 1.25;                      % our target for aggregate markup
        
        N            = 143.24;                    % mean number of firms in sector
        p.sigma      = 1.1485;
        p.gamma      = 7.1573; 
        p.xi         = 5.1451; 

    case '1.35'    

        Markuptarget = 1.35;                        % our target for aggregate markup
         
        N            = 112.44;                      % mean number of firms in sector
        p.sigma      = 0.99424;
        p.gamma      = 5.2127; 
        p.xi         = 3.7245; 

end

% Assigned Parameters

p.beta       = 0.96;                        % period is 1 year
p.varphi     = 0.04;                        % exit rate 
p.delta      = 0.06;                        % depreciation rate

p.r          = 1/p.beta - 1; 
p.nu         = 1;                           % inverse labor supply elasticity
p.alpha      = 1/3;                         % capital share in value added
p.theta      = 0.5;                         % elasticity of substitution between VA and intermediates
 
S            = 25000*4;                     % number of sectors

Nmax         = max(100, floor(N*2));        % maximum number of firms in sector

% Pareto tail productivity

Y            = 1; 
R            = p.r + p.delta; 

mshare       = 0.45;                        % share of intermediate inputs in total sales

% Quadrature
 
nz          = 100;                          % Gaussian discretization for z

zmax        = -1/p.xi*log(1e-12);

[p.z, p.w]  = qnwlege(nz, 0, zmax); 

p.w         = p.xi.*p.w.*exp(-p.xi.*p.z);

p.w         = p.w./sum(p.w); 
p.z         = exp(p.z); 

cut         = cumsum(p.w) > 1 - 1e-5; 

wend        = sum(p.w(cut)); 
zend        = p.w(cut)'*p.z(cut)/wend; 

p.w         = [p.w(~cut); wend]; p.w = p.w/sum(p.w);
p.z         = [p.z(~cut); zend]; 

nz          = size(p.z, 1); 
Fcum        = [0; cumsum(p.w)];              % cumulative distribution

Fcum(end)   = 1; 

p.z         = p.z/3; 

% Random draws

rng(0); 

% 1. Allocate firms to sectors

ns          = max(1, poissrnd(N, S, 1)); 
Nf          = Nmax*S; 

% 2. Productivity draw from discretized distribution

unif        = nodeunif(Nf, 0, 1); 
unif        = unif(randperm(Nf)); 

[~, ~, bin] = histcounts(unif, Fcum);

z           = reshape(p.z(bin), Nmax, S); 
zsave       = z; 
nssave      = ns; 

Nsave       = mean(ns); 

for i = 1 : S
   
    z(ns(i) + 1 : end, i) = eps;         % kill extra firms (set z = eps)
    
end

nsmax         = max(ns); 
z             = z(1 : nsmax, :);         % no reason to carry extra entries

z             = sort(z, 1, 'descend');

% Solve equilibrium markups in every sector

[mu, prof, omega, Mus, Zs, Omegas, Mu, Z, Omega]  = industryequilibrium([], z, p);

p.phi            = 1 - mshare*Omega^(1 - p.theta)*Mu;                                       % choose to match materials share in sales = mshare 

W                = (1 - p.alpha)*(((Omega^(1 - p.theta) - (1 - p.phi))/p.phi)^(1/(1 - p.theta))*(R/p.alpha)^(-p.alpha))^(1/(1 - p.alpha));        % Invert the CES Cost aggregator                                              % Labor only factor of production                

Lp               = (1 - p.alpha)*p.phi/Mu*((R/p.alpha)^p.alpha*(W/(1 - p.alpha))^(1 - p.alpha)/Omega)^(1 - p.theta)*Y/W; 
K                =       p.alpha*p.phi/Mu*((R/p.alpha)^p.alpha*(W/(1 - p.alpha))^(1 - p.alpha)/Omega)^(1 - p.theta)*Y/R;
B                = (1 - p.phi)/Mu*(1/Omega)^(1 - p.theta)*Y;

C                = Y - p.delta*K - B;

moment_data        = [Markuptarget; 0.45; 0.430; 0.723; 0.213]; 

moment_model       = zeros(size(moment_data)); 

moment_model(1)    = Mu; 
moment_model(2)    = B/Y;

sales              = omega.*Omegas*Y;
labor              = (1 - p.alpha)*p.phi*((R/p.alpha)^p.alpha*(W/(1 - p.alpha))^(1 - p.alpha)/Omega)^(1 - p.theta)*sales./mu;  % W*l

lshare             = sum(labor)'./sum(sales)';               % sectoral labor share 
hh                 = sum(omega.^2)'; 

bhat               = lscov([ones(S, 1), hh], lshare);

CR4                = Omegas*sum(omega(1:4, :))'/sum(Omegas);
CR20               = Omegas*sum(omega(1:20, :))'/sum(Omegas);

moment_model(3)    = CR4;
moment_model(4)    = CR20;
moment_model(5)    = -bhat(2);

weights        = ones(numel(moment_model), 1); 

weights(1)     = 100; 
weights(5)     = 10; 

weights        = weights/sum(weights);

err_mom        = (moment_model - moment_data)./(1 + moment_data);
err_mom        = (weights'*err_mom.^2).^(1/2);

price   = mu./z*Omega;
Ps      = sum(price.^(1 - p.gamma)).^(1/(1 - p.gamma)); 
P       = mean(Ps.^(1 - p.sigma))^(1/(1 - p.sigma)); 

Ys      = (Ps./P).^(-p.sigma).*Y;
y       = (price./Ps).^(-p.gamma).*Ys; 

% Efficient Allocations

Zse     = sum(z.^(p.gamma - 1)).^(1/(p.gamma - 1)); 
Ze      = mean(Zse.^(p.sigma - 1)).^(1/(p.sigma - 1));

Musave  = Mu; 
Zsave   = Z; 
Zesave  = Ze; 

GDPe  = p.phi^(1/(p.theta - 1))*Ze/(1 - (1 - p.phi)*Ze^(p.theta - 1))^(1/(p.theta - 1))*K^p.alpha*Lp^(1 - p.alpha);
Ye    = p.phi^(1/(p.theta - 1))*Ze/(1 - (1 - p.phi)*Ze^(p.theta - 1))^(p.theta/(p.theta - 1))*K^p.alpha*Lp^(1 - p.alpha);
Be    = (1 - p.phi)*Ze.^(p.theta - 1)*Ye; 

GDP   = p.phi^(1/(p.theta - 1))*(1 - (1 - p.phi)*Z^(p.theta - 1)*Mu^(-p.theta))*Z/(1 - (1 - p.phi)*Z^(p.theta - 1)*Mu^(1 - p.theta))^(p.theta/(p.theta - 1))*K^p.alpha*Lp^(1 - p.alpha);
GDP1  = p.phi^(1/(p.theta - 1))*(1 - (1 - p.phi)*Z^(p.theta - 1)              )*Z/(1 - (1 - p.phi)*Z^(p.theta - 1)                 )^(p.theta/(p.theta - 1))*K^p.alpha*Lp^(1 - p.alpha);

weight   = price.*y./mu;   % proportional to each firm's labor

% Sectoral Stats

weights  = sum(labor)';              % sectoral labor shares

% Output

clc

switch report

    case 'T5'

        fprintf('<strong> Table 5: Parameterization, Oligopoly </strong> \n'); 

        fprintf('\n');
        
        fprintf([' Aggregate Markup at ', num2str(round(Mu,2)), '\n'])

        fprintf('\n');
        
        symbol      = categorical({'M'; 'CR4';'CR20';'-';'b';sprintf('%c ', 958);sprintf('%c ', 947);sprintf('%c ', 951);'N';sprintf('%c ', 952)});
        description = categorical({'aggregate markup';'top 4 sales share';'top 20 sales share';'materials share';'regression coefficient';'Pareto tail';'elasticity of substitution within sectors';'elasticity of substitution between sectors';'average number of firms per sector';'weight on valued-added'});
        values      = [Mu; moment_model(3); moment_model(4); B; -moment_model(5); p.xi; p.gamma; p.sigma; round(N,0); p.phi];
        
        disp(table(symbol, description, round(values, 2), 'VariableNames', {'Symbol','Variable Description','Value'}))

    case 'T6'

        fprintf('<strong> Table 6: Markup Dispersion and Productivity Losses, Oligopoly </strong> \n'); 

        fprintf('\n');
        
        fprintf([' Aggregate Markup at ', num2str(round(Mu,2)), '\n'])

        fprintf('\n');
        
        description = categorical({'25th percentile';'50th percentile';'75th percentile';'90th percentile';'99th percentile'});
        values1      = [wprctile(Mus(:), 25, weights); wprctile(Mus(:), 50, weights); wprctile(Mus(:), 75, weights); wprctile(Mus(:), 90, weights); wprctile(Mus(:), 99, weights)];
        values2      = [wprctile(mu(z > .1), 25, weight(z > .1)); wprctile(mu(z > .1), 50, weight(z > .1)); wprctile(mu(z > .1), 75, weight(z > .1)); wprctile(mu(z > .1), 90, weight(z > .1)); wprctile(mu(z > .1), 99, weight(z > .1))];
        
        fprintf(' Cost-weighted distribution of markups \n');  
        
        fprintf('\n');
        
        disp(table(description, round(values1, 2), round(values2, 2), 'VariableNames', {'Variable Description','Sector Level','Firm Level'}))
        
        fprintf('\n');
        
        description = categorical({'gross output ';'value added';'value added, M = 1'});
        values      = [-(Z/Ze - 1)*100; -(GDP/GDPe  - 1)*100; -(GDP1/GDPe - 1)*100];
        
        fprintf(' Aggregate productivity loss in percent \n');  
        
        fprintf('\n');
        
        disp(table(description, round(values, 2), 'VariableNames', {'Variable Description','Value'}))
        
        fprintf('\n');

    case 'T7'

        % Compute expected profits by integrating pi(z, s) 
        
        Profz       = zeros(nz, 1); 
        
        Znew        = zeros(nz, 1); 
        Munew       = zeros(nz, 1); 
        
        parfor i = 1  : nz
            
           muold    = [repmat(p.gamma/(p.gamma - 1), 1, S); mu]; 
        
           znew     = [repmat(p.z(i), 1, S); z]; 
            
           [~, profnew, ~, ~, ~, ~, Munew(i), Znew(i)]  = industryequilibrium(muold, znew, p);
        
           Profz(i) = mean(profnew(1, :)); 
        
        end
        
        p.F        = (p.w'*Profz)*Y/W/(1/p.beta - 1 + p.varphi);
        
        G          = (p.w'*Profz);
        Q          = Y*G/(1 - p.beta*(1 - p.varphi));
        
        % Find total employment
        
        L          = p.F*p.varphi*N + Lp; 
        
        p.psi      = W/L^p.nu/C; 

        clc

        switch experiment

            case {'uniform', 'entry'}

                start_subsidy

            case {'efficient', 'sizedep'}

                start_planner

        end

end
